Spectral function of the Anderson impurity model at finite temperatures 
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Using the functional renormalization group (FRG) and the numerical renormalization group 
(NRG), we calculate the spectral function of the Anderson impurity model at zero and finite temper- 
atures. In our FRG scheme spin fluctuations are treated non-perturbatively via a suitable Hubbard- 
Stratonovich field, but vertex corrections are neglected. A comparison with our highly accurate 
NRG results shows that this FRG scheme gives a quantitatively good description of the spectral 
line-shape at zero and finite temperatures both in the weak and strong coupling regimes, although at 
zero temperature the FRG is not able to reproduce the known exponential narrowing of the Kondo 
resonance at strong coupling. 
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I. INTRODUCTION 

In the past decade, the experimental realization of 
man-made nanostructures such as quantum dots coupled 
to a metallic environment^ has brought renewed atten- 
tion to the theoretical handling of the Anderson impurity 
model (AIM). Moreover, the solution of the AIM repre- 
sents one of the fundamental steps in the so-called dy- 
namical mean-field theory^ describing the quantum dy- 
namics of the Hubbard model in the limit of infinite spa- 
tial dimensions^ where the effects of the dynamical Weiss 
field can be described in terms of an effective local im- 
purity model subject to a self-consistency condition. Be- 
cause the thermodynamics of the AIM can be obtained 
exactly by means of the Bethe Ansata^ and, on the 
other hand, Wilson's numerical renormalization groupl 
(NRG) gives a numerically controlled method of calculat- 
ing spectral properties^ - — the AIM can also be used as 
a benchmark model for testing various non-pcrturbative 
many-body methods. It is therefore important to de- 
velop reliable methods for solving the AIM especially at 
intermediate and strong coupling which require low com- 
putational effort, in the perspective of tackling the more 
complex impurity problems arising from both dynamical 
mean-field theory applications and experimental realiza- 
tions of nanostructure devices. 

Motivated by the above considerations, several authors 
have tried to reproduce at least some aspects of the 
known properties of the AIM using functional renormal- 
ization group (FRG) methods ii2r-i£ Although the FRG 
gives a formally exact hierarchy of integro-differential 
equations for all irreducible vertices^ - — in practice this 
hierarchy has to be truncated in order to extract physical 
information from it. Unfortunately, at this point no trun- 
cation of this hierarchy has been found which correctly re- 
produces all the known strong-coupling properties of the 
AIM, such as the correct interaction dependence of the 
Kondo scale, which is known to be exponentially small 
in the interaction strength. In fact, the proper construc- 
tion of unbiased truncations of the formally exact FRG 
flow equations which do not break down in the strong 



coupling limit is a largely unsolved problem in the field. 

Quite generally, some progress in solving the FRG flow 
equations can be made if the dominant fluctuation chan- 
nel is known a priori. For example, for the description 
of superfluidity in the attractive Fermi gas the particle- 
particle channel is known to play a special role. In such 
a situation it is natural to decouple the interaction in the 
particle-particle channel using a Hubbard- Stratonovich 
transformation and consider the FRG flow of the cou- 
pled mixed fermion-boson modeh 19 ' 21- — 

In Ref. [l5| a similar strategy has been applied to 
the AIM, although in this case the proper choice of the 
Hubbard-Stratonovich decoupling is not so obvious. For 
simplicity, we shall focus in this work on the local mo- 
ment regime of the symmetric AIM at intermediate to 
strong coupling; in this case the physics is dominated by 
spin fluctuations, so that one should decouple the on- 
site interaction in the spin-fluctuation channels. Still, 
the decoupling is ambiguous because of the freedom of 
distributing the interaction between transverse and lon- 
gitudinal spin channels. In the simplest case, the decou- 
pling is then performed such that only transverse spin- 
fluctuations are explicitly introduced. The resulting FRG 
flow equations have already been derived in Ref. [l5j . 
but further approximations have been made to reduce 
the resulting integro-differential equation to an ordinary 
differential equation for the wave-function renormaliza- 
tion factor Z . Unfortunately, in this approximation only 
the low-energy features of the spectral function can be 
described. In this work we shall not rely on such a low- 
energy approximation, but present a fully self-consistent 
solution of the integro-differential equation for the self- 
energy of the AIM. This enables us to calculate the full 
spectral line-shape of the AIM at all energy scales. In 
order to test the accuracy of our FRG approach, we shall 
also calculate the spectral function at zero and finite tem- 
peratures using the highly reliable NRG approach. By 
comparing the results obtained with the two methods, we 
then show that our rather simple truncation of the FRG 
flow equations gives a quantitatively accurate description 
of the spectral function of the AIM at all energies. 
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II. FRG WITH PARTIAL BOSONIZATION IN 
THE SPIN SINGLET CHANNEL 

The AIM describes a single localized electron level 
which is coupled to a band of non-interacting conduc- 
tion electrons. The latter degrees of freedom can be inte- 
grated out so that the partition function can be written 
as a fcrmionic functional integral involving only Grass- 
mann fields d a (r) and d a (r) associated with the localized 
electrons. Here a =t, 4- labels the spin projection, and 
r is the imaginary time. Following Rcf. [l5| we decou- 
ple the local on-site interaction by means of a complex 
bosonic Hubbard-Stratonovich field x i n the spin-singlet 
particle-hole channel. The partition function Z divided 
by the corresponding partition function -Ehf m Hartree- 
Fock approximation can then be written as 



Zuf fT>[d,d,x,x]e~ So[d ' d 'X' x] 
The Gaussian part of the bare action is 



(1) 



So[d,d,x,x] = ~ y^G 1 (iu)d u<T d ua - 
Ju a 

+ f U^x^Xu,, (2) 
while the interaction part can be written as 

Si[d,d,x,x] = / [du+ufduiXu + dwidu+utXu] , 

where J = (3 1 ^ w denotes summation over fermionic 
Matsubara frequencies iuj, with (3 = 1/T being the in- 
verse temperature, whereas J Q = Z? -1 ^^ denotes sum- 
mation over bosonic Matsubara frequencies iuj (in the 
zero temperature limit the sums are replaced by inte- 
grals). For simplicity we consider only the symmetric 
AIM in the wide band limit, where the Hartree-Fock 
propagator is given by 



G (iui) 



1 



iA sign a; 



(4) 



In the above expression the energy scale A arises from 
the hybridization between the localized d-clcctrons, rep- 
resented by the Grassmann variables d CT and d a , and the 
conduction electrons. The exact d-electron propagator, 
on the other hand, is of the form 



G(iw) 



iA sign a; — 



(5) 



Our aim is to calculate the irreducible self-energy T,(iuj) 
using the FRG and then analytically continue it to the 
real frequency axis. 

We use here the frequency transfer cutoff scheme pro- 
posed in Ref. [TEj , where an infrared cutoff is introduced 
only into the bosonic part of the Gaussian action ©. 



Specifically, we introduce the cutoff via the following sub- 
stitution in the second line of Eq. ^ , 



where 



ir 



R A {iuj) 



U- 1 + R A (iu), 



A 

7rA 



R(\u,\/A), 



and the function R(x) is defined by2£ 

R{x) = (l-x)O(l-x). 



(0 



(7) 



(8) 



Our implementation of the FRG method is therefore dif- 
ferent from previous purely fermionic FRG studies of the 
AIM; 13 i 14 where no Hubbard-Stratonovich fields were in- 
troduced. The advantages of our scheme are that the 
resulting FRG flow equations can be analyzed with mod- 
erate numerical effort and yield a quantitatively accurate 
description of both the low-energy and high-energy fea- 
tures (such as the broadened Hubbard peaks) of the spec- 
tral function. Neglecting vertex corrections, the FRG 
flow of the cutoff-dependent self-energy reads 

<9a£ a («w) = / F A (iul)GA(iuj - iuj), (9) 



where the flowing fermionic propagator G A (iui) depends 
again on the flowing self-energy, 



1 



(10) 



iuo + iA sign a; — Y, A (iuj) ' 
and the bosonic single-scale propagator is given by 

F A fa) = [-d A R A (iLo)][F A (iu>)} 2 , (11) 
with the bosonic propagator given by 

1 



F A (iQ) 



U- 1 + R A (iuj) ~U A (iQ)' 



(12) 



Here, H A (iu~j) represents the flowing irreducible spin sus- 
ceptibility, which can be related to the fermionic Green 
function, in the absence of vertex corrections, by means 
of the following (skeleton) equation^ 



(13) 



H A (iuj) = — G A {iuj)G A (iuj — iuj) 



Substituting Eqs. <jTU]) - CC2) into Eq. © we obtain a 
closed integro-diffcrential equation for the flowing self- 
energy ~E A (iuj), which must be integrated from the initial 
scale Ao — > oo, with initial condition T, Ao (iuj) = 0, down 
to A = 0. The physical self-energy is then given by the 
limit T,(iuj) = limA->o ^a(*w). 

Before discussing the details of our numerical approach 
to Eq. ©, it is worthwhile to remark that in Ref. [ll| 
the corresponding flow equation for Y. A (iuj) was treated 
in a very crude approximation, namely by expanding all 
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quantities appearing in Eq. ((9j) to linear order in fre- 
quency. As a consequence, only the low-energy proper- 
ties of the spectrum, characterized by the ui — quasi- 
particle peak, were accessible. Moreover, the results were 
obtained only in the zero-temperature limit. Here, on the 
other hand, we do not rely on any low-energy expansion, 
so that we can calculate the spectral function at all en- 
ergy scales and describe the transfer of spectral weight, 
in the strong-coupling regime, from the Kondo peak to 
the high-energy Hubbard bands. 



III. NUMERICAL SOLUTION OF THE FRG 
EQUATION FOR THE SELF-ENERGY 

A. Zero temperature 

In the zero temperature limit all the sums over Mat- 
subara frequencies are replaced by integrals, /3 _1 ^ w — > 
J j^. The flow equation for Sa(«w) reduces therefore to 
an integro-diffcrcntial equation, which requires an arti- 
ficial discretization of the imaginary axis in order to be 
solved numerically. More precisely, our numerical ap- 
proach to Eq. ([5]) consists in replacing the flowing self- 
energy H\(iuj), defined as a function of the continuous 
parameter u), by a finite set of flowing couplings Ea(«w,i) 
defined on a discrete mesh of frequencies w„. 

The choice of a specific paramctrization of the Mat- 
subara axis, at zero temperature, is in principle arbitrary; 
however, one has to carefully verify that the numerical re- 
sults are actually independent of such a parametrization. 
In other words, the results should be numerically stable 
with respect to the choice of the discretization procedure. 
To achieve this goal we consider the following geometric 
meshi^ of frequencies, 

a n - 1 

0J n = Wmin — , U = 1,...,N, (14) 

a — 1 

where the free parameters a > 1 and oj m i n > define the 
spacing of the frequencies, while N is the total number of 
frequencies to be used (note that in the particle-hole sym- 
metric case we have H\(iu)) = — Sa(— ioj), so that we do 
not need independent couplings for the negative frequen- 
cies). As thoroughly discussed in Ref. [14} . the unequal 
spacing of such a mesh allows indeed to resolve with great 
accuracy the low-energy regime of the spectrum, which 
is known to exhibit the sharp Kondo resonance at uj = 0, 
while covering, at the same time, a sufficiently wide range 
of frequencies in order to rule out the spurious effects of 
a finite frequency cutoff w max = ^>n- More specifically, 
one has to simultaneously satisfy, for all the values of U 
under consideration, the following conditions, 



w miD < Z{U)A, 
Wmax > max(A, U), 



where 



Z(U) = 1 - 



d(iuj) 



(15) 
(16) 



(17) 



is the quasiparticlc residue (proportional to the width of 
the Kondo peak), and afterwards verify that the results 
are numerically stable upon a further increase in both 
the spanned frequency range [w m in,w max ] and the num- 
ber N of sampling frequencies. A typical choice of the 
discretization parameters is, e.g., 



= l(T b A, 



1.06, N = 400, 



(18) 



which allows to solve the flow equations with a relatively 
small computational effort up to the strong-coupling 
regime U < 87rA. 

The resulting system of N coupled differential equa- 
tions for the flowing couplings Y,^(iuj n ) is solved by means 
of refined Runge-Kutta routines, using a linear interpo- 
lation procedure in order to evaluate the Green function 
in Eqs. ([H]) and (|13p at frequencies that do not belong to 
the geometric mesh (i.e., when cj£ — Gift ^ {u„}). 



B. Finite temperatures 

At finite temperatures the Matsubara space is intrin- 
sically discrete, with fermionic and bosonic frequencies 
given by 



oj n = (2n + l)7r//3, u n = 2mr/P, 



(19) 



so that the flowing self-energy itself is already defined as 
a countable set of couplings T,A(iu! n ). Eq. © is therefore 
no longer a proper integro-diffcrcntial equation, as in the 
zero temperature limit, but consists rather of an infinite 
set of coupled ordinary differential equations, which can 
be solved numerically by keeping only the first N Mat- 
subara frequencies. 

In contrast to the zero temperature case, it is crucial to 
observe that now the spacing between consecutive points 
of the mesh is constant, its value 6u> = 2irT being dic- 
tated by the actual value of the temperature T. This fact 
simplifies the implementation of the numerical routines, 
since there is no need for any interpolating procedure; on 
the other hand, however, it limits the possibility of study- 
ing arbitrarily small temperatures. In fact, in order to 
obtain numerically stable results using the natural Mat- 
subara discretization defined in Eq. (|T9"|) . one typically 
needs N > 10 3 for temperatures T/A < 0.1, which al- 
ready corresponds to a remarkable computational effort. 
In other words, although smaller values of T/A are still 
accessible from the numerical point of view, the computa- 
tional time required becomes comparable to that of fully 
numerical methods (e.g. the NRG), making the present 
method numerically unworthy. 



C. Analytical continuation 

The solution of our flow equation gives, by construc- 
tion, the physical self-energy along the imaginary axis, 
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H(iu)) = £a=o In order to access real frequency 

properties, such as the spectral function 

A(u) = --ImG(iw->w + iO), (20) 

we need to analytically continue our results for cither 
G(iuj) or E(iw) onto the real frequency axis. Since our 
numerical solutions are not affected by statistical errors, 
we perform the analytical continuation using the Pade 
approximation,— which we found to give more stable 
results if directly applied to the self-energy. In other 
words, we first evaluate the real frequency self-energy, 
= Yi{iuj — > lu + iO), and afterwards calculate the 
corresponding Green function using the Dyson equation. 
The frequency dependence of the self-energy is indeed 
much more sensitive to the interaction-induced features 
of the model than the Green function itself, whose non- 
trivial behavior is masked by the asymptotic form of the 
non-interacting propagator. 



IV. NUMERICAL RENORMALIZATION 
GROUP AT FINITE TEMPERATURES 

The numerical renormalization group method^ was 
first applied to the Anderson impurity model^ shortly 
after its introduction in 1975. Originally designed to 
calculate thermodynamic quantities at very low temper- 
atures, the NRG has been extended over the last two 
decades to calculate more complex observables, such as 
the single particle spectral function investigated here (for 
a recent review covering many extensions and applica- 
tions of the NRG, see Ref. [l2|]). We briefly mention 
several crucial developments which lead to improving the 
accuracy of dynamical correlation functions calculated by 
NRG: the z-trick averaging^ which improves the resolu- 
tion at frequencies close to the conduction-band edge; the 
introduction of the reduced density matrbs^ in order to 
account for the correct low-tempcrature state when cal- 
culating the spectrum at higher frequencies; and finally, 
a recently developed smart choice of the basis for the full 
NRG Fock spaced which was designed to calculate the 
real-time evolution after a quantum quench, and succes- 
sively has been applied to improve the accuracy of spec- 
tral functions for the single impurity Anderson models 
by avoiding overcounting of spectral transitions within 
the NRG. 

The NRG data presented here were obtained by the 
full density matrix approach^ which gives access to re- 
liable finite-temperature spectra for frequencies as low as 
uj ~ T/5. In combination with the self-energy trick 31 
and an average over N z = 32 slightly different discretiza- 
tions of the conduction band^ we obtain high-quality 
NRG spectra, which serve as a benchmark for the newly 
developed FRG scheme discussed above. A detailed de- 
scription of the NRG method is beyond the scope of this 
publication; for details we refer the interested reader to 
Ref. nlf and references therein. 
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FIG. 1: (Color online) Spectral function of the symmetric 
AIM evaluated with the FRG (solid lines) and the NRG 
(dashed lines) for different interactions and temperatures. 
From top to bottom the temperature increases from T/A = 
to T/A = 0.2, as indicated. 



V. RESULTS 

Our results for the spectral function at temperatures 
T/A = 0, 0.05 and 0.2 arc plotted in Fig. [Q where in 
each panel we have considered several values of the in- 
teraction strength, ranging from the weak to the strong 
coupling regime. The zero temperature line-shapes are 
in good agreement with the NRG data in the weak and 
strong coupling regimes, while in the intermediate cou- 
pling range U / (7rA) m 1 our method somewhat overesti- 
mates the role of the interaction, enhancing the transfer 
of spectral weight from the Kondo resonance to the Hub- 
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FIG. 2: (Color online) FRG spectral function for different in- 
teractions and temperatures. From top to bottom the interac- 
tion increases from U/(ttA) = 1 to U/(nA) = 6, as indicated. 
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FIG. 3: (Color online) Inverse quasiparticle weight Z _1 at 
zero temperature on a logarithmic scale, as a function of 
U/(nA). The inset shows the weak-coupling behavior in the 
regime < U/(tvA) < 1 on a larger scale. The present FRG 
method corresponds to the solid curve, the dashed and dot- 
dashed curves are the FRG results of Ref. [l5l |. the black cir- 
cles represent the NRG data and the dot-dot-dashed curve in 
the inset is the second order perturbation theory. 



consideration. This behavior is more clearly shown in 
Fig. [2j where the FRG results are now plotted, in each 
panel, for a fixed value of the interaction and increasing 
temperatures. 

Finally, in Fig. [3] we show the inverse quasiparticle 
weight at T = as a function of the interaction. Unfortu- 
nately, in contrast to the NRG results, at strong coupling 
our FRG method is not able to reproduce the exponen- 
tial narrowing of the quasiparticle peak predicted by the 
Bethe Ansatz£ However, the numerical estimate of \jZ 
for large values of U / (ttA) is improved compared to the 
previous FRG results presented in Ref. [l5j]. Fig. [3] shows 
also clearly the above mentioned overestimation of the 
interaction in the intermediate coupling regime, charac- 
terized by a stronger renormalization of the quasiparticle 
peak. 



bard bands. It is worthwhile to observe that for large 
values of the interaction the peaks of the Hubbard bands 
are located approximately at ±[//2 and their width is 
of order 2A, in good agreement with the known strong- 
coupling limit results^ As the temperature increases, the 
Kondo resonance peak becomes more and more broad- 
ened, as expected when the temperature approaches the 
Kondo scale Tk- However, we notice that at small and 
intermediate couplings this effect is more enhanced in the 
FRG approach, in comparison to the NRG data. On the 
other hand, the position and the shape of the Hubbard 
bands remain essentially unchanged, being related to an 
energy scale much larger than the temperatures under 



VI. CONCLUSIONS 

In the present work we have calculated the spectral 
function of the symmetric Anderson impurity model at 
zero and finite temperatures using functional and nu- 
merical renormalization group methods. In particular, 
we took advantage of our high quality NRG data in or- 
der to test the validity of the FRG truncation scheme 
proposed in Ref. [15j . where the on-site interaction is 
decoupled in the spin-singlet particle-hole channel by 
means of a bosonic Hubbard-Stratonovich field and an in- 
frared cutoff is explicitly introduced only into the bosonic 
propagator. Following Ref. [l5[, we neglected the FRG 
flow of vertex functions in order to truncate the in- 
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finite hierarchy of FRG flow equations and obtain a 
closed integro-differential equation for the flowing self- 
energy £a(ju;). In contrast to Ref. [IH, however, where 
only the low-energy expansion of the self-energy was 
taken into account, we solved numerically the integro- 
differential equation for Ea(^) keeping the complete fre- 
quency structure. This has allowed us to calculate the 
spectral function at all energy scales, including the high- 
energy Hubbard bands. 

Comparing the FRG results to the NRG data, we 
found that our truncation scheme gives quantitatively 
good results for the spectral line-shape in the weak and 
strong coupling regimes, both at zero and finite temper- 
atures, although it somewhat overestimates the effects of 
the interaction at intermediate couplings [//(ttA) rs 1. 
Most important, in contrast to purely fermionic FRG 
approachc o 13 ! 14 where no Hubbard-Stratonovich fields 
are introduced, our method gives an accurate descrip- 
tion of both the low-energy and high-energy features of 
the spectral function, including the correct position and 
width of the high-energy Hubbard bands in the strong 
coupling regime. 



Unfortunately, our FRG truncation scheme is still not 
able to correctly reproduce, at zero temperature, the ex- 
ponential narrowing of the Kondo peak in the strong cou- 
pling regime. Most likely, such a non-perturbative effect 
requires exact symmetry relations, e.g. Ward identities^ 
to be enforced throughout the integration of the FRG 
flow. 
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